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ABSTRACT 

The Fermi Gamma-ray Space Telescope reveals two large bubbles in the Galaxy, which extend 
nearly symmetrically ~ 50° above and below the Galactic center (GC). Using three-dimensional (3D) 
magnetohydrodynamic (MHD) simulations that self-consistently include the dynamical interaction 
between cosmic rays (CR) and thermal gas, and anisotropic CR diffusion along the magnetic field lines, 
we show that the key characteristics of the observed gamma-ray bubbles and the spatially-correlated 
X-ray features in ROSAT 1.5 keV map can be successfully reproduced by a recent jet activity from the 
central active galactic nucleus (AGN). We find that after taking into account the projection of the 3D 
bubbles onto the sky, the physical heights of the bubbles can be much smaller than previously thought, 
greatly reducing the formation time of the bubbles to about a Myr. The 'young' bubbles naturally 
satisfy the upper limit of bubble ages estimated from the cooling time of high-energy electrons. No 
additional physical mechanisms are required to suppress large-scale hydrodynamic instabilities because 
the evolution time is too short for them to develop. The simulated CR bubbles are edge-brightened, 
which is consistent with the observed projected flat surface brightness distribution. Furthermore, we 
demonstrate that the sharp edges of the observed bubbles can be due to anisotropic CR diffusion along 
magnetic field lines that drape around the bubbles during their supersonic expansion, with suppressed 
perpendicular diffusion across the bubble surface. Possible causes of the slight bends of the Fermi 
bubbles to the west are also discussed. 



1. INTRODUCTION 

One of the most important findings from the first 
two years of Fermi Gamma-ray Space Telescope obser- 
vations are two large bubbles extending to ~ 50° above 
and below the Gal actic center (GC), with a width o f 
~ 40° in longitude (jSu et al.ll2010l; IDobler et al.|[20Toh . 
The Fermi bubbles are nearly symmetric about the GC, 
with only slight bends toward the west (negative longi- 
tude). They have approximately flat gamma-ray sur- 
face brightness with sharp edges. The bubbles emit 
at 1 < < 100 GeV, and have a spatially-uniform 
hard spectrum {dN 1 /dE 1 ~ E~ 2 ). The gamma-ray 
emission is most likely to originate from inverse Comp- 
ton (IC) scattering of photons in the interstellar radi- 
ation field (ISRF) by cosmic ra y (CR) electron s with 
energies 10 < E cr < 100 GeV (jSu et all l2010h . The 
same population of CR electrons is invoked to ex- 
plain the hard-spectrum microwave synchrotron radi- 
ation spatially-correlated with the Fermi bubbles ob- 
served by the Wilkinson Microwave Anisot ropy Probe 
(WMAP), also known as the 'WMAP haz e' (fFinkbcincr 
120041: IDobler fe Finkbeinerl 120081 : iDoblerl I2012T ). The 
edges of the bubbles also coincide with features in the 
ROSAT X-ray maps at 1.5 keV (jSnowden et all Il997t 
ISu et alJ[2010h . 

As discussed in ISu et al.l (|2010l) , the unique location 
and morphology of the Fermi bubbles suggest that they 
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are created by some large episode of energy injection from 
the GC, such as a nuclear starburst in the last ~ 10 Myr, 
or a past accretion event onto the central supermassive 
black hole (SMBH). The latter scenario, i.e., the forma- 
tion of CR-filled bubbles by recent jet activity of the 
active galactic nucleus (AGN), is appealing due to sev- 
eral reasons. Relativistic jets from AGN can accelerate 
CR electrons to high energies and thus directly provide 
the source for the gamma-ray emission. The hard, spa- 
tially uniform spectrum of the Fermi bubbles requires 
the CR electrons to be transported from the sites where 
they are generated to where the bubbles are observed 
today without significant cooling, which means the bub- 
bles must be at most a few million years old (i.e., the 
IC cooling time of ~ 100 GeV electrons). Assuming 
the cosmic rays are produced at the GC, to travel to 
a distance of several kpc within a few Myr, they must 
be transported very rapidly, at a speed of ^transport ~ 
10 4 (Z/10 kpc)(i ago /l Myr)- 1 km s~\ which is readily 
achievable by the fast AGN jets. Furthermore, AGN 
jets are known to be responsible for radio synchrotron 
emiss ion from ext e nded extragalactic rad io sources 
(e.g., IScheuerl Il97l IBlandford fc Reesl Il97l and CR- 
filled bubbles (e.g.. lLaing et al.ll2006fl observed in mas- 
sive g alaxies and galaxy clusters ( McNamara fe Nulsenl 
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120071) . though their gamma-ray emission cannot be easily 
detected due to limited sensitivity and resolution. The 
Fermi bubbles could possibly be an analogy to such cases, 
and their proximity may provide a special opportunity to 
study AGN bubbles in the gamma-ray band. 

Forming the Fermi bubbles in this jet scenario has re- 
cently b een explored using two-d imensional (2D) simula- 
tions bv lGuo fc Mathews! (|20 111) . They found that AGN 
jets are able to efficiently carry the cosmic rays to sev- 
eral kpc within ^1 — 3 Myr, and the axial ratios of 
the bubbles can be reproduced when the density con- 
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trast of the jets relative to the surroundings is within 
the range 0.001 < 77 < 0.1. Despite the general success, 
there are several discrepancies with the observed bub- 
bles. Their simulated bubbles are subject to large-scale 
hydrodynamic instabilities that induce ripples on the side 
of the bubbles, in contrast to the rather smooth surface 
of the observed bubbles. Also, their simulated bubbles 
have a uniform CR distribution, which would appear to 
be limb-darkened if projected onto the sky. In order to 
ameliorate these problems, they proposed shear viscosity 
as a possible mechanism to suppress the instabilities as 
well as to produce an edge-brightened CR distribution 
that is necessary for a flat projected surface brightness 
profile (|Guo et al.N20Tl . 

They also showed that the sharpness of bubble edges 
requires suppression of CR diffusion across the bubble 
surface; otherwise, the bubble edges would be substan- 
tially smoothed if the cosmic rays dif fuse isotropically 
with typical coefficient s in the Galaxy ()Guo fc Mathews! 
1201 It IGuo et al.ll201lh . The key to account for the re- 
quired suppression of CR diffusion may be the effect of 
anisotropic diffusion, i.e., diffusion of cosmic rays along 
magnetic field lines with strongly inhibited cross-field dif- 
fusion, expected when the gyro-radius of the cosmic rays 
is much smaller (~ 10~ 9 kpc for a 10 GeV electron in a 4 
/iG magnetic field) than the mean free path between colli- 
sions. If during the expansion of the CR bubbles the am- 
bient gas is compressed, resulting in tangential magnetic 
field lines near the surface, then the suppression of dif- 
fusion across the bubble boundaries could be explained. 
Using hydrodynamic simul ations with specified sp atial 
variation in CR diffusivity, IGuo fe Mathews! (|201 lh and 
IGuo et all (|2011f ) showed that the observed sharp edges 
can indeed be produced if the diffusion coefficient at the 
bubble surface is significantly smaller than the canoni- 
cal values. However, as also stressed by these authors, 
understanding how this process occurs requires detailed 
modeling of both the magnetic field and anisotropic CR 
diffusion along the field lines. 

To this end, we perform three-dimensional (3D) mag- 
netohydrodynamic (MHD) simulations of the formation 
of the Fermi bubbles by CR jet injections from the 
central SMBH in the Galaxy. Our simulations self- 
consistently include the effects of magnetic field, CR ad- 
vection, dynamical coupling between cosmic rays and 
thermal gas, and CR diffusion along the field lines. The 
main objective of this study is to use 3D numerical sim- 
ulations involving additional physical mechanisms to in- 
vestigate whether the jet scenario is able to produce bub- 
bles that are consistent with the observed features, in- 
cluding the shape, the flat surface brightness, and the 
sharp edges. We will highlight the c ompa risons with the 
previous work of IGuo fc Mathews! ll20ll regarding the 
differences between pure hydrodynamics and MHD, and 
between 2D and 3D simulations, which have significant 
impact on the predicted bubble morphology. Moreover, 
we will show that the interplay between anisotropic CR 
diffusion and magnetic fields can account for the sharp 
edges of the observed Fermi bubbles. The structure of 
this paper is as follows. In §[2] we describe the numerical 
techniques and initial conditions employed. In § [31 we 
first present characteristics of our simulated CR bubbles 
and compare them to the morphology of the observed 
bubbles, and then study the influence of magnetic fields 



and anisotropic CR diffusion on the sharpness of bubble 
edges in § 13.21 We note the parameter dep endencies and 
available observational constraints in § 14.11 discuss possi- 
ble causes of the slight bends of the observed bubbles in 
§ 14.21 and compare our results with X-ray observations 
in § 14.31 Finally, the summary and implications of our 
findings are given in § [5j 

2. METHODOLOGY 
2.1. Assumptions and Numerical Techniques 

We perform 3D MHD simulations of CR injections 
from the GC usi ng the adaptive-mesh-refinement (AMR ) 
code FLASH4 (|Frvxell et alJUM IDubev et all 12008ft . 
The simulation box is 50 kpc on a side, refined progres- 
sively on the GC located at the center of the simula- 
tion domain. The resulting coarsest and finest resolution 
elements are 0.8 and 0.1 kpc, respectively. The diode 
boundary condition is used, which is similar to the out- 
flow boundary condition but does not allow matter to 
flow into the domain. The MHD equations are computed 
using t he directionally unsplit stagg ered mesh (USM) 
solver (|Lee fc Deand 120091 : Eeel 120121 ). The USM algo- 
rithm in FLASH4 is based on a finite- volume, high-order 
Godunov scheme and is combined with a constrained 
transport method to ensure the divergence-free condi- 
tion of magnetic fields. The order of the USM algorithm 
used in our simulations corresponds to the P i ecewis e- 
Parabolic Method fPPM; ICoIella fc Wood ward! (fl98l ). 
which is well suited for capturing shocks and contact 
discontinuities, as needed in our simulations. In order 
to incorporate cosmic rays in the simulations, we imple- 
mented a new CR module in the USM solver in FLASH4. 
Here we give a brief summary of the assumptions and 
equations used to simulate cosmic rays, and present the 
implementation details and numerical tests of the CR 
module in Appendix [A"l 

Though each CR particle travels near the speed of 
light, it is well known that their collective transport 
speed in the Galaxy is much smaller, due to scatter- 
ing by magnetic irregularities. When the scattering is 
significant, cosmic rays are trapped by the magnetic ir- 
regularities that are frozen in the plasma. The cosmic 
rays effectively move with the thermal gas at the local 
gas velocity, i.e., are 'advected'. Cosmic rays may stream 
through thermal gas in regions where the magnetic field 
is locally alig ned, with streaming speed lim ited by the 
Alfven speed (jSkillind Il97lt IKulsrudl 12001 . Since the 
Alfven speed is much smaller than the local gas speed 
(as the jet is supersonic), in this study we will neglect 
the effect of CR streaming. Cosmic rays can also diffuse 
with respect to the thermal gas as they scatter off mag- 
netic inhomogeneities. Though CR diffusion with a typ- 
ical diffusivity in the Galaxy is too slow to be the main 
transpor t mechanism for the cosmic rays in th e Fermi 
bubbles (ISu et al.ll2010l : IGuo fc Mathews! [20T1 . it can 
have non-negligible effects on the CR distribution within 
and at the edges of the bubbles and thus is included in 
our simulations. 

To simulate collisionless cosmic rays together with col- 
lisional plasma, a complete description would require 
solving for the distribution function of cosmic rays, 
f(r,p,t). Because / depends on three momentum co- 
ordinates as well as position and time, the system of 



3 



equations is numerically more difficult than a set of fluid 
equations. In many situations where / is nearly isotropic 
in momentu m space, / can be solved as a functio n of only 
r, and t l|Skilung||1975t iMiniati et al.ll200l . In this 
paper we adopt a more simplistic approach, which is to 
treat cosmic rays as a second fluid and solve directly for 
th e evolution of CR pressure p cr as a function of r and 
t (IDrurv fe Voelkl[T98lt IJones fe K an gl [19901 IRvu et al 



2003; Mathews fc Brighcnti 2008; Rascra fe Chandran 



2008). The cosmic rays are advected with the thermal 
gas, and in return the gas can react to the gradients of 
the CR pressure. In this approach, the cosmic rays are 
treated as a single species without distinction between 
electrons and protons, though CR electrons may be the 
primary source of the gamma-ray emission of the Fermi 
bubbles. We did not model the CR energy spectrum, 
and neglected the cooling and heating processes of cos- 
mic rays, such as energy losses due to synchrotron and 
IC emission, and reacceleration in shocks. We will inves- 
tigate their effects in a future work. 

The MHD equations including CR advection, diffusion, 
and dynamical coupling between the thermal gas and 
cosmic rays can be written as 



dp 
dt 
dpv 



V • (pv) = 0, 



pvv 



BB 

47T 



Vptot = pg, 



3B 

— -Vx(vxB)=0, 



de 
dt 



(e + Ptot)v - 



B(Bv) 



4tt 



(1) 
(2) 
(3) 

(4) 



= pv ■ g + V • (k ■ Ve cr ), 
<9e cr 

— ^- + V • [e CI v) = -p cr V • v + V • (« • Ve cr ), (5) 
at 

where p is the gas density, v is the velocity, B is the mag- 
netic field, g is the gravitational field, n is the diffusion 
tensor, e cr is the CR energy density, and e = 0.5pv 2 + 
eth + e cr + B 2 /8tt is the total energy density. The total 
pressure is p to t = (7 - l) e th + (7cr - l)e cr + B 2 /8tt, where 
eth is the internal energy density of the gas, 7 = 5/3 is 
the adiabatic index for ideal gas, and 7 cr = 4/3 is the 
effectiv e adiabatic index of co smic rays in the relativistic 
regime (|Jubelgas et al.ll2008[ ). 

In the presence of a magnetic field, the diffusion term 
can be written as 



V • (k ■ Ve cr ) = V • (n\\bb ■ Ve cr ) + 
V-K(J-66)-Ve cr ], 



(G) 



where km is the diffusion coefficient parallel to the mag- 
netic field, k± is the perpendicular diffusion coefficient, 
and b is the unit vector of the magnetic field ( Braginskii 
[19651) . 

In discussions of Galactic CR propagation, km and k± 
are defined with the respect to the mean magnetic field, 
not the exact field (see also § 12. 2p . Perpendicular diffu- 
sion is assumed to be due to some combination of field- 
line wandering due to a random field component and 
magnetic field structure on scales less than the CR gy- 
roradius, which allows the cosmic rays to migrate from 



one fieldline to another (iDuffv et al.lll995t UokipiilH99l 
lEnBlin fe Voetl [20031 IHauff et all I2010D . In our imple- 
mentation, we recognize that the magnetic field is prob- 
ably structured on scales below our numerical resolution, 
so kii and k± are necessarily defined with respect to 
a mean. In the simulations presented here, we assume 
that the true field is oriented along the mean field, and 
so we set k± = 0. Since in typical cases the perpendic- 
ular diffusion coefficient is m uch smaller than the par- 
allel diffusion coefficient (e.g., Enfilin fc Vogt 2003]), our 
test results show that including a nonzero but small per- 
pendicular diffusion coefficient does not change our main 
conclusions. 

When magnetic field lines are sufficiently tangled on 
small scales, CR diffusion can be approximated to be 
isotropic, for which the diffusion term in the equations 
can be expressed as Ki so V 2 e cr . Assuming isotropic dif- 
fusion, typical values of the diffusion coefficient in the 
Galaxy is observationally constrained to be Kj SO ~ (3 — 
5) x 10 28 cm 2 s , which ma y depend on the m agnetic 
field structure and CR energy (jStrong et al.H2007j ). In or- 
der to compare the effects of isotropic and anisotropic CR 
diffusion, especially on the sharpness of the Fermi bub- 
bles, we perform simulations for both cases with varied 
diffusion coefficients. Since in a tangled field the effective 
diffusion coefficient (Ki so ) can be suppressed compared to 
that along the field lines (km), likely by a factor of ~ 3 
(|Taoi ri995). we choose to vary the parallel diffusion co- 
efficient in the range of (0.4 — 1.2) x 10 29 cm 2 s _1 . The 
simulation parameters are summarized in Table [TJ 

2.2. The Galactic Model 

In order to form the Fermi bubbles, we inject CR 
jets from the GC into pre-existing gas in the Galactic 
halo. X-ray observations of emission or absorption 
lines have provided evidences for the existence of 
hot gas with temperature T ~ 10 6 K extending out 
to the virial radius of the Milky Way halo with a 
scale height of a few tens of kpc dBlitz fc Robishawl 
20001; iMcCammon et"all 120021; iRasmussen et all 1 2003 1 ; 



Fang et alj 12 006 
Grcevich fc Putmanl 



Bregman fc Llovd-Daviea 



2007 



.2009 1 ; IMiller fc Bregman] 120121) . 
Following IGuo fc Mathewsl (|2011D . we initialize the hot 
gaseous halo to be in hy drostatic equilibrium within a 
fixed Galactic potential (|Helmi fc White! I2001D . assum- 
ing the gas is initially isothermal. In our 3D simulation 
domain, the Galactic potential can be written as 



* = $halo + $disk + $bulgc, 



where 



$halo = fhalo + d\) 

is the dark logarithmic halo, 

GAfdisc 



^x 2 + y 2 + (a + \J1FT~B 1 ) 2 

is the Miyamoto-Nagai disk, and 

GA/buigc 
"Pbulgc — — ; — 



(7) 
(8) 

(9) 
(10) 



is the spherical Hernquist stellar bulge. The distance 
to the GC is r = \J x 2 + y 2 + z 2 , and the GC is in the 
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TABLE 1 
Simulation Parameters 



Run Magnetic Field ?b (kpc) CR Diffusion 

^iso OI" ^f| 

(cm 2 s i) 



A 


None 




None 




B 


Tangled 


9 


None 




C 


Tangled 


9 


Isotropic 


4.0 X 10 28 


D 


Tangled 


9 


Anisotropic 


4.0 x 10 28 


E 


Tangled 


9 


Isotropic 


1.2 x 10 29 


F 


Tangled 


9 


Anisotropic 


1.2 x 10 29 


G 


Tangled 


1 


None 




H 


Tangled 


1 


Isotropic 


4.0 x 10 28 


1 


Tangled 


1 


Anisotropic 


4.0 x 10 28 


J 


Tangled 


1 


Isotropic 


1.2 x 10 29 


K 


Tangled 


1 


Anisotropic 


1.2 x 10 29 



middle of the computational domain extending from —25 
kpc to +25 kpc in x-, y-, and z-directions. The adopted 
parameters for the potential are fhaio = 131.5 km s -1 , 
d h = 12 kpc, M disc = 10 n M Q , a = 6.5 kpc, b = 0.26 
kpc, Afbuigc = 3.4 x 10 10 M Q and d h = 0.7 kpc. 

The gas density is related to the electron number den- 
sity by 

p = p,en e m^, (11) 

where m M is the atomic mass unit, /i e = 5/x/(2 + fi) is 
the molecular weight per electron, and fi — 0.61 is the 
molecular weight. Given the electron number density 
at the origin, n e o, and the initial gas temperature, To, 
the gas density distribution is derived from hydrostatic 
equilibrium. We choose n e o = 2 cm~ 3 and Tq = 2 X 10 6 
K in order to match the o bserved hot gas density profile 
([Miller fc Bregmanl ([201 2D ; see also § Q)> . We neglect 
rotation in the galactic disk because the rotation period 
is much longer compared to the estimated age of the 
Fermi bubbles. 

The hydrostatic model presented here is not the only 
possible one. The hot gas distributio n in the inner galaxy 
is als o well fit by a galactic wind ([Everett et al.l I2008L 
I2010T ) emanating from a ring with inner galactocentric 
radius about 3 kpc. Although such a wind would have 
no effect on the initial expansion of the bubbles, it would 
provide a rather different background medium in the 
later stages. We plan to explore this in future work. 

One of the main objectives of this paper is to investi- 
gate anisotropic CR diffusion caused by magnetic fields in 
our Galaxy. To this end we initialized the ambient fields 
(as opposed to injected magnetic fields from the jets, see 
§ 12. 3[) based on available observational constraints. The 
magnetic field in our Galaxy is composed of a large-scale 
regular field (or the mean field) and a small-scale tur- 
bulent field. The latter is associated with the turbulent 
interstellar medium, which has a typical coherence leng th 
of ~ 5 — 50 pc (see a recent review by iNouts'o s (2012)). 
Since this scale is not resolved in our simulation, we will 
only model explicitly the regular field but not the tur- 
bulent component. Therefore, we note that CR diffusion 
considered in our simulations should be interpreted as an 
effective diffusion with respect to the mean field. 

Observational constraints on the strength and struc- 
ture of the global Galactic magnetic field (GMF) have 
been greatly improved over the past decade; however, 
some details are still unclear. For example, the magnetic 
fields in the galactic disk is estimated to be a few /u,G and 



enhanced in the spiral arms and near the GC. The direc- 
tion of the disk field roughly follow directions of the arms, 
but the numb er of field reversals is still under debate (see 
iBrownl (|2010( ) and references therein). The halo field also 
has a typical magnitude of ~ 1 fiG. Nevertheless, there 
has been no agreement on a dipole, quadrupole, o r even 
a large-scale coherent vertical field (Brown 2010). Due 
to the uncertainties in the detailed magnetic-field con- 
figurations, instead of imposing a certain shape for the 
regular field, we model the GMF with a tangled field 
with random orientations. In this paper we show results 
from two sets of simulations with different magnetic field 
coherence lengths, i.e., 9 kpc (as in the halo) and 1 kpc 
(as in the disk), in order to bracket possible cases where 
the coherence length is larger or smaller than the size of 
the Fermi bubbles, respectively. 

The initial tangled field is computed outside the main 
simulation code and is generated by 3D inverse Fourier 
transform (FFT) of ma gnetic field that, in fc-sp ace, has 
an amplitude given by ([Ruszkowsk i et al.ll2007 " 



B k oc fc- 11/6 exp(-(fc/fc 1 ) 4 )exp(-fc 2 /fc), 



(12) 



where k = (k%. + k 2 + k 2 ) 1 / 2 , and k\ and A)2 are the cutoff 
wavenumbers for large and small k, respectively. Given 
the 3D magnetic power spectrum, B{k) oc |-Bfc| 2 , the co- 



herence length can be defined as I b 
coher 
2003) 



2tt Iks, where the 

coher ent scale in fc-space is defined as (Enfilin fc Vogtl 



f 



k 2 B{k)dk 



/„ kB(k)dk 



(13) 



In our simulations with Is = 9 kpc, we set k\ = 0.92 and 
&2 = 0.3; for the case of Ib — 1 kpc, k\ = 6.98 and k% = 3 
are used. To ensure V • B = 0, the following projection 
is performed for 'cleaning' the field of divergence terms 
in Fourier space: 



B'(k) = (I-kk)B(k), 



(14) 



where J is the identity operator, and k is the unit vec- 
tor in fc-space. We note that this operation does not 
change the power spectrum of the magnetic field fluctu- 
ations. After performing the inverse FFT, the magnetic 
field strength is normalized such that the average field 
strength is 1 fiG. The initial magnetic field for these two 
cases in the central region is shown in the left panels of 
Figure [3l 
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2.3. Jet Injection 

The feasibility of producing the Fermi bubbles by AGN 
jets and parame t er var iations are studied in detail by 
iGuo fe Mathews! (|201lD . In order to facilitate compar- 
isons, and in particular to focus on the effects of mag- 
netic fields and anisotropic CR diffusion, we set up our 
jet model (for thermal gas and CR components) based 
on their prescriptions but include a few modifications. 
We also inject magnetic fi elds along with the jets using 
the AGN subgrid model of lSutter etall (|2012t) . O ur nu- 
merical scheme of jet injection is as follows. 

The bipolar jets are introduced at the GC along the 
+z and —z direction with a constant velocity Vj e t for a 
duration of £j e t from the start of the simulation. The 
direction of the jets are parallel to the rotational axis of 
the Galaxy, motivated by the symmetry of the observed 
northern and southern bubbles (except the case in § 14.21 
where we suggest an explanation for the slight bends of 
the observed bubbles) . The thermal gas and CR contents 
of the jets can be described using six parameters: the 
thermal gas density p 3 , the gas energy density ej , the CR 
energy density ej cr , the jet velocity i>j et , the radius of the 
jet cross-section i?j C t, and the jet duration tj Ct . The sum 
of kinetic, thermal, and CR power is thus 

^jet = -Pkc + ^th + Per = (^Pivf et + ej + e jcr ^ nRf et Vj et . 

(15) 

In IGuo fc Mathews! (|201 If ) . the jets are imposed in a 
cylinder inside the simulation domain. However, this 
method would introduce an inconsistency between the 
distance jets can travel in one timestep and the height 
of the cylinder, and an artificial redistribution of jet en- 
ergy inside the cylinder. To avoid these drawbacks, we 
inject the jets through a disk on the z = plane with 
radius i?j e t using the inflow boundary condition, i.e., up- 
dating the ghost cells (in the middle of the computa- 
tional domain) with the jet parameters at each simu- 
lation timestep before solving the MHD equations. In 
this way variables in the simulation domain are updated 
self-consistently according to the exact fluxes across the 
boundary calculated from the USM solver. 

Jets from SMBH are generally thought to carry mag- 
netic field, as it plays an important role in jet collima- 
tion and acceleration fsee iFerra ri (1998) and references 
therein), and as suggested by observations of lobes of 
radio galaxies (lOwen et al.l 120001 : IKronberg et al.l [20011 : 
ICrostonet al.H2005T) . AGN are also considered to be a vi- 
able agent for distributing the fields into the intergalactic 
medium (IGM) an d intracluster medium (ICM) by rel- 
ativistic jets (e.g., IKoide et al.l 119991 : IContopoulos et all 
120091 : ISutter et al.l 12012^ A realistic representation of 
the GMF thus requires modeling of the magnetized jets. 
The details of launching the jets in the relativistic regime 
is beyond the scope of our study; here we only consider 
the large-scale effects of the magnetized jets after they 
have propagated to the resolvable scale. To this end, we 
inject magneti c fields using the subgrid magnetized jet 
model of ISutter et al.1 (|2012D . This model is based on 
the magnetic 'tower' model of iLi et all ([20061), in which 
the magnetic field is composed of a poloidal flux that pre- 
sumably comes from the dynamo process in black hole 
accretion disks, as well as a toroidal component that may 



be generated by shearing of the poloidal flux lines by dif- 
ferential rotation in the disks. T he spatial configura tion 
of the injected field is described in lSutter et al.1 (|2012t ). to 
which we refer the readers for more details. During the 
active phase of the jet, a fraction Jb of Pj' et is injected 
in the form of magnetic energy onto the grid as a source 
term in the induction equation (right hand side of Eq. 
[3]). The /b is essentially a free parameter, and there- 
fore in our simulations it is set such that the resulting 
magnetic field in the end of the simulations is consistent 
with observed values. We find t hat the observ ed field 
strength of > 10 pG at the GC (|Noutsosl [2012]) can be 
reproduced using fg = 10~ 3 , which means the jets are 
not magnetically dominated. 

A few additional quantities are used to characterize 
the jets: the density contrast between the thermal gas 
contained in the jets and the ambient gas, r\ = Pj/Vamb, 
the thermal energy contrast ry = ej/e am b, the total jet 
power Pjot = (1 + /s)Pj' ct , and the total injected energy 
Ej ct = Pj e tijet- Note that since we inject the jets from 
the GC, the ambient gas density and internal energy are 
defined at the initial values at t he origin, as opposed to a 
few kpc away from the GC as in IGuo fe Mathe^i(|20Tl . 
Due to the cuspy gas density profiles near the GC (see 
their Figure 1), our derived p am b and e am b are higher 
than their quoted values. Therefore, the jets in our sim- 
ulations in general are more powerful than theirs, so as 
to overcome the ambient gas pressure near the origin. 

We have verified that whe n using identical i mplem en- 
tations and parameters as in IGuo fc Mathewi(|201lT) . wo 
are able to reproduce their results. However, as also 
pointed out by the authors, the observed shape of the 
Fermi bubbles can be recovered by multiple parameter 
combinations due to degeneracies. Moreover, the results 
(e.g., the derived jet power) are dependent on the ini- 
tial gas profiles employed. Therefore, in this study we 
make an effort to incorporate available observational con- 
straints on the jet parameters and initial gas profiles 
(see § 12.21 and § 14. ip . By observing X-ray absorption 
lines through the hot g aseous halo along many d ifferent 
sight lines in the sky, iMiller fc Bregmanl (|2012f) found 
that the ratio between OVIII and OVII column densities 
is enhanced for sight lines that pass through the Fermi 
bubbles, indicating gas temperature of a few times 10 
K inside the bubbles, much lower than the temperature 
produ ced by the fiducial parameters in IGuo fc Mathewi 
(|2011l ). We find that matching the observed bubble tem- 
perature requires slower and wider jets, together with 
adjusted density contrast and CR energy density in or- 
der to maintain the observed shape of the bubbles. To 
this end, we choose a set of fiducial jet parameters for 
the simulations presented in this paper. Their values are 
given in the upper part of Table [2j Given the fiducial pa- 
rameters and the initial profile of the Galactic halo, the 
characteristic quantities of the jets are derived and listed 
in the bottom part of Table [2J In § 14.11 we will discuss 
the parameter dependence of bubble morphology, and 
we will show that there are in fact very limited degrees 
of freedom for the jet parameters given all the available 
observational constraints on the properties of the Fermi 
bubbles. 

3. RESULTS 
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TABLE 2 
Input and Derived Jet Parameters 



Parameter 


Description 


Value 


Unit 


V 


Density contrast 


0.05 




Ve 


Energy density contrast 


1.0 




e jcr 


CR energy density 


2.5 X 10 -9 


erg cm " 


''jet 


Jet speed 


0.025 


c 


Rjet 


Radius of cross-section 


0.5 


kpc 


tjet 


Duration of injection 


0.3 


Myr 


n ej 


Electron number density 


0.1 


cm ° 


Pi 


Thermal gas densiy 


1.95 x 10" 25 


g cm -3 


e i 


Thermal energy density 


1.59 x 10~ 9 


erg cm - 3 


1'he 


Kinetic power 


3.08 x 10 44 


erg s 


Pth 


Thermal power 


8.90 x 10 42 


erg s 


Per 


CR power 


1.40 x 10 43 


erg s — 1 


1'b 


Magnetic power 


3.31 x 10 41 


erg s — 1 


Piet 


Total power 


3.31 x 10 44 


erg s — 1 




Total injected energy 


3.13 x 10 57 


erg 



The total injected energy by both bipolar jets is 2E 



'jet- 



3.1. Morphology of Fermi Bubbles 

We first present simulations of the Fermi bubbles with- 
out magnetic field and CR diffusion (i.e., Run A), fo- 
cusing on the compar isons with the previous work of 
IGuo fc Mathews! (|2011l ). namely, the differences between 
2D and 3D simulations. In particular, we find that pro- 
jections of the 3D bubbles are nontrivial and have signif- 
icant impacts on the interpretations of bubble properties 
previously derived. 

Firstly, the estimation of the formation time of the 
Fermi bubbles are significantly influenced by the consid- 
eration of the projection effect. In the simulations, the 
formation time of the bubbles can be defined as the time 
when the vertical extent of the simulated bubbles reaches 
the observed height of |6| ~ 50°. Assuming the bub- 
ble sizes are negligibl e compared to the distan ce between 
the Sun and the GC, IGuo fc Mathews! (|2011l ) found that 
this condition is equivalent to a vertical bubble size of 
~ 10 kpc in the simulation domain, and the estimated 
age of the bubbles is generally around 1 — 3 Myr. In our 
simulations, the coordinate transformation from the 3D 
simulation domain (x, y, z) to the 2D projected map in 
Galactic coordinates (l,b) is 



tan^ : 



tan b = z 



V + Rq 

y + Ro 



cos I 



(16) 
(17) 



where Rq = 8 kpc is the distance from the Sun to the 
GC, and the bubbles are observed from the Sun's coor- 
dinate, defined to be (x,y,z) — (0, — i? Q ,0). If we were 
to adopt the time when the bubbles reaches \z\ = 10 
kpc as the bubble age, the projected map of the bub- 
bles would have an extent of \b\ ~ 60° — 70°, because 
the widths of the bubbles (~ 4 — 5 kpc) are actually 
comparable to Rq, and therefore the front side (closer 
to the observer) of the 3D bubbles would be projected 
to a higher galactic latitude and reach the observed lati- 
tude much earlier than what would be naively expected. 
Taking into account the effect of projection, we find that 
the condition of \b\ ~ 50° corresponds to \z\ ~ 6 kpc. 



Therefore, we stop the simulations when the height of 
the CR bubbles of e cr > 10~ n erg cm~ 3 along the jet 
axis reaches 6 kpc, and define the time t = t Fermi as the 
current age of the Fermi bubbles. Due to the projection 
effect, the age of bubbl es we derive is much sho rter than 
previous estimation by IGuo fc Mathews! (|2011D . For in- 
stance, when we use the same initial conditions and jet 
parameters as theirs, the age of the bubbles is only ~ 1/3 
of their estimation, solely due to the consideration of 
projection (note that since the bubble expansion veloc- 
ity slows down over time, it takes much shorter time for 
the bubbles to expand from kpc to 6 kpc than from 6 
kpc to 10 kpc). 

The short age of the bubbles derived has several im- 
portant implications. Firstly, the Fermi bubbles are ob- 
served in the energy range of 1 < E 1 < 100 GcV, which 
corresponds to relativistic CR electrons with 10 < E CI < 
100 GeV if the gamma-ray emission is produc ed by iC 
scatt ering of the ISRF by these CR electrons (|Su et al.1 
2010). Since the IC cooling time of CR electrons at the 
high-ener gy end [~ 100 GeV) is only a few Myr (Fig- 
ure 28 in ISu et "all (|2010D ). this puts an upper limit on 
the age of the Fermi bubbles. In our AGN jet-inflated- 
bubble scenario with the consideration of projection, it 
only takes 1.2 Myr for the cosmic rays to travel 6 kpc 
from the GC to the observed latitude (note that our jets 
are slower than those used in lGuo &F Mathews 120111 and 
the bubble formation time would have been much longer 
if the bubbles had to reach 10 kpc). Therefore, the con- 
straint from the IC cooling time is naturally satisfied by 
our 'young' bubbles. Later in the section, we will show 
that the short bubble formation time also has critical in- 
fluence on the dynamics and morphology of the bubbles 
as observed today. 

Figure [1] shows the CR energy density and the ther- 
mal electron number density sliced through the y = 
plane for the fiducial run A at tFermi = 1.2 Myr. In this 
pure-hydrodynamic run, the bubbles are symmetric with 
respect to the z = plane, and thus only the north- 
ern bubble is plotted. Initially the jet is over-pressurized 
with respect to the ambient medium and energetically 
dominated by the kinetic power. The cosmic rays in- 
jected by the central AGN form a pair of CR bubbles 
above and below the GC. In this run CR diffusion is not 
included, and thus the evolution of the bubbles is solely 
due to CR advection and the effects of CR pressure. Due 
to the high velocity of the thermal gas in the jets, the 
cosmic rays are rapidly advected to z = 6 kpc at the end 
of the simulation iFcrmi = 1.2 Myr. The lateral expansion 
of the bubbles, on the other hand, is mainly driven by the 
large pressure contrast between the bubbles (contributed 
by both the thermal gas and the cosmic rays) and the 
ambient pressure declining with the distance from the 
GC. Note that the dent on the bubble surface along the 
jet axis is a result of the cuspy initial gas density pro- 
file. Because there is a non-negligible density gradient 
across the width of the jets, it is more difficult for the 
central part of the jets to penetrate the dense ambient 
material at the GC. If a more core- flattened initial den- 
sity profile were chosen, the top of the bubbles would be 
flatter and reach 6 kpc a little earlier, causing a slightly 
shorter iFcrmi- However, because the lateral expansion is 
relatively unaffected, the shape of the projected CR bub- 
bles, which is mainly determined by the side edges closer 
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Fig. 1. — Slices of CR energy density (top) and thermal electron 
number density (bottom) in simulation coordinates for the fiducial 
run A at <F cr mi = 1.2 Myr. Only the northern bubble is plotted 
since the southern bubble is symmetric with respect to z = 0. 
Quantities are shown in logarithmic scale and in cgs units. 



to the Sun, remains unchanged, and only the projected 
CR distribution within the bubbles is slightly affected. 

As can be seen from the bottom panel of Figure [TJ 
the hot halo gas is expelled and compressed into a shell 
due to shocks produced by the fast expansion of the CR 
bubbles. Inside the shocks the electron number density 
is enhanced to n G ~ 0.02 cm -3 and the gas is heated 
to T < 10 8 K, separated by a contact discontinuity with 
the hot, underdense gas with n c = 10~ 4 ~ 5 x 10 -3 cm~ 3 
and T > 10 7 K inside the bubbles. At i Fermi = 1.2 Myr, 
the shocks are moving at a speed of 6000 — 7000 km s _1 
(Mach number M ~ 30) at the shock front along the z 
direction and ~ 2500 km s _1 (Mach number M ~ 12) 
in the lateral direction at z = 3 kpc. Due to com- 
pression and heating, the shocked gas has an enhanced 
bremsstrahlung emissivity and result in limb-brightened 
X-ray emission around the bubbles (see § 14.31) . The 
strong shocks are in contrast to weak shocks (M ~ 
1 — 2) associated with radio bubbles and X-ray cavities in 
galaxy clusters. Moreover, while in clusters buoyancy is 
often important in the evolution of radio bubbles, it plays 
a lesser role for the simulated Fermi bubbles as their ex- 
pansion is dominated by the momentum and energy in- 
jected by the jets. Note that our prediction of shock ve- 
locities at the prese nt day is a few times larg er than pre- 
vious estimation by IGuo fc Mathewsl (|2011h . Again be- 
cause of the short bubble formation time, the shocks are 
observed today at an earlier epoch of the evolution and 
have not slowed down as much in the ambient medium. 
The time for the shock front to move one arcsecond (re- 
solvable with the Chandra X-ray Observatory) in the lat- 
eral direction is ~ 13 yr. 
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Fig. 2. — Projected maps of CR energy density in Galactic co- 
ordinates for the fiducial run A at tp e «ni = 1.2 Myr, plotted in 
logarithmic scale in cgs units. The solid and dotted lines show 
the surfaces of the observed northern and southern Fermi bubbles, 
respectively. 

The projected maps of CR energy density in Galactic 
coordinates for the fiducial run A at i Fermi = 1.2 Myr are 
shown in Figure The intensity on this map is roughly 
proportional to the gamma-ray surface brightness assum- 
ing the gamma-ray emission is produced by IC scattering 
by CR electrons of a smoothly distributed ISRF. As can 
be seen from the plot, the projected CR distribution after 
evolving for i Fermi = 1.2 Myr resembles the shape of the 
observed Fermi bubbles (solid and dotted line s for the 
northe rn and southern bubbles, respectively; ISu et al.l 
( 2010)), with only slight discrepancies on the east side 
because of the bends of the observed bubbles (see § 14.21 
for the discussion and modeling of bubble bending). One 
of the important characteristics of the observed bubbles 
is their sharp edges. We find that in the absence of CR 
diffusion (as in Run A) , the sharpness is also seen for the 
simulated gamma-ray bubbles, implying suppression of 
CR diffusion across the bubble surfaces. In § 13.21 we will 
show that the suppression can possibly be attributed to 
CR diffusion along magnetic field lines that drape around 
the bubble surface. 

It is nontrivial to reproduce the flat gamma-ray sur- 
face brightness of the observed Fermi bubbles. Assum- 
ing the ISRF photon density is roughly constant or de- 
creases away from the GC, the CR energy density must 
increase toward the bubble boundaries in order to yield 
a uniform distribution after line-of-sight projections. As 
shown in the top panel of Figure [TJ our simulated CR 
energy density indeed is smaller in the bubble interior 
and gradually increases toward the rims, which then is 
projected to produce a rather uniform surface bright- 
ness (see Figure [2|) , especially along the lateral direction 
for a fixed Galactic latitude. The simulated CR energy 
density after projection does have a dependence on the 
Galactic latitude b (greater at a higher latitude). Elec- 
tron aging is an unlikely explanation for the observed 
flat brightness distribution as a function of latitude, be- 
cause within the short bubble formation time, IC cooling 
of CR electrons of energy 10 — 100 GeV is negligible at 
> 1 kpc a way f rom the Galactic plane (see Figure 28 of 
ISu et al.l (|2010h ). Instead, since the ISRF is expected 
to be stronger near the Galactic plane, the theoretical 



IC radiation distribution, which is proportional to the 
product of the CR energy density and ISRF, may be 
flatter as a function of b than the emission naively in- 
ferred just from the projected CR distribution as a func- 
tion of b. Furthermore, the observed surface brightness 
distribution is subject to larger errors closer to the Galac- 
tic disk due to background subtraction. Therefore, the 
observed flat distribution of gamma-ray emission as a 
function of latitude is not inconsistent with our simula- 
tions. We note that th e fiducial parameters adopted in 
IGuo k, Mathews! fj201 lh produce a uniform simulated e cr 
distribution, which implies a limb- darkened gamma-ray 
surface brightness after projection, contrary to the flat- 
ness of the observed surface brightness. For this reason, 
these authors suggested physical viscosity as a possible 
mechanism to produce a CR distribution cons istent with 
the gamma-ray observation (|Guo et al.1 120111 ). Here we 
show that by using jet parameters appropriate for an 
observationally-constrained initial gas density profile (see 
§ l4.1l for discussion on the chosen parameters), it is possi- 
ble to reproduce the observed uniform surface brightness 
without introducing additional physical mechanisms. 

Another feature to note is the rather smooth surface 
of the observed Fermi bubbles, which is also reproduced 
for our fi ducial jet paramet e rs. Th is is in contrast to the 
results of IGuo fc Mathewsl (|2011l ). who show that their 
simulated bubbles are rippled due to large-scale hydro- 
dynamic instabilities. We find that the difference may 
lie in the estimation of bubble formation time. As dis- 
cussed in the beginning of this section, our bubble ages 
are generally shorter than their estimations because it 
takes much less time for the projected 3D bubbles to 
match the observed bubble sizes. Given the young age 
of the bubbles, the hydrodynamic instabilities, including 
Rayleigh- Taylor (RT) and Kelvin-Hehnholtz (KH) insta- 
bilities, do not have enough time to develop before the 
bubbles are observed today (though the instabilities do 
show if the bubbles are observed at later times). When 
only the effect of gravity is considered, the timescale for 
the growth of RT instability at the bubble surface can be 
estimated by 

t _ 12 , Myr ^ ( ^ r (_^)''\ 

(18) 

where ?7 S is the density contrast across the contact dis- 
continuity at the bubble surface, g s is the magnitude 
of gravitational acceleration in cgs units, and A s is the 
wavelength of the mode of instability under considera- 
tion. For modes comparable to the radius of the bub- 
bles (A s ~ 3 kpc) at the bubble surface on the jet axis 
(g s (x = 0,2/ = 0,z = 6 kpc) - 1(T 8 , rj B ~ 0.25), the 
RT timescale is £rt ~ 15.8 Myr. This estimate is likely 
only a lower limit, because the condition for the onset of 
RT i nstability becomes more stringent for inflating bub- 
bles (|Pizzolato fc Sokerl 120061) . Firstly, the deceleration 
of bubbles would counteract the gravitational accelera- 
tion in the frame of reference comoving with the bubble 
surface in the early stages of the evolution. Also, the RT 
growth rate would slow down even further because of an 
additional drag force as well as an effect of perturbation 
stretching associated with bubble inflation. Therefore, 
the RT instability grow on a timescale much longer than 



the ages of the bubbles, and hence it has a negligible 
effect. 

On the other hand, the KH timescale is more relevant 
since it is comparable to the estimated ages of the bub- 
bles 

(19) 

where the shear velocity Av s and density contrast rj s are 
scaled to typical values at the bubble surface along the 
lateral direction (e.g., x = 3 kpc, y = 0, z = 3 kpc). If 
£ Fermi > ^kh, the KH instability would grow efficiently 
and produce ripples in the CR distribution, as found in 
IGuo &: Mathewsl (|201 If) : otherwise there is no sufficient 
time for the instability to develop, and the resulting bub- 
bles should have a smooth surface, as seen in our simu- 
lations. Since the ripples are in apparent contradiction 
to the smooth surface of the observed bubbles, these au- 
thors propos e to add viscosity as a means to suppress the 
instabilities ()Guo et al.l 1201 If ). By taking into account 
the effect of 3D projection, we show that the formation 
time of the Fermi bubbles are in fact shorter than previ- 
ously expected, which can naturally explain the absence 
of hydrodynamic instabilities and the smooth surface of 
the observed bubbles, even when magnetic fields are ne- 
glected. 

3.2. Effects of CR Diffusion on Bubble Edges 

One of the most important characteristics of the ob- 
served Fermi bubbles is the sharp edges of their sur- 
face brightness. If CR diffusion is purely isotropic with 
a diffusion coefficien t Ki so ~ (3 — 5) x 10 28 cm 2 s _1 
(| Strong et al.l 120071 ). the resulting CR distribution 
would have a much smoother edge than observed 
(|Guo & Mathewsl 120111 ). The sharpn ess of bubble edges 
thus imp ly suppression of CR diffusion across the bubble 
surfaces (|Guo fc Mathewsll201lHGuo et al.ll2011f) . Using 
MHD simulations that self-consistently include effects of 
magnetic fields on CR diffusion, we will show that this 
suppression is likely due to anisotropic CR diffusion along 
magnetic field lines that drape around the Fermi bubbles. 

Figure [3] shows the magnetic field strength and topol- 
ogy at t = and t = tFermi for runs with different initial 
field coherence length, Run B (Zb = 9 kpc; top row) and 
Run G (Ib = 1 kpc; second row). A zoom-in image of 
Run G at iFermi is shown in the bottom panel for clarity. 
For both simulations, the initial tangled field has an av- 
erage strength of ~ 1 fiG, which corresponds to plasma 
ft = Ptot / (B 2 / 8ir) ~ 1 in the ambient galactic halo (note 
that p cr = for the ambient gas), in agreem ent with ob- 
servations (e.g.. IHaverkorn fc Heesenll2012D . During the 
supersonic expansion of the Fermi bubbles, the ambi- 
ent gas is compressed into shells with enhanced density 
and m agnetic fields. Because of the magnetic draping 
effect (|Lvutikovl 120061 ). the field lines are stretched and 
compressed in the draping layer, causing field alignment 
with the bubble surface as well as amplification of field 
strength. The resulting alignment and amplification are 
greater when the direction of the field is initially parallel 
to the bubble surface (perpendicular to the direction of 
propagation). At t = i Fermi, the magnetic field within the 
shells is amplified to > 10 /iG, comparable to the equipar- 
tition field strength inferred by radio polarization mea- 
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Fig. 3. — Slices in simulation coordinates of magnetic field magnitude (in logarithmic scale in units of /^G) and directions (shown in 
arrows) at £ = (left) and t = tF orm i (right) for simulations without CR diffusion, Run B (top row) and Run G (second row). Run B and 
G has an initial magnetic coherence length of 9 kpc and 1 kpc, respectively. For clarity, a zoom-in image of Run G at iFermi is shown in 
the bottom panel. Due to the effect of magnetic draping, the field is amplified and aligned with the bubble surface. 



10 



surements (|Jones et al.ll2012T) . Depending on the initial 
field strength and direction, the degree of field alignment 
and amplification at t — t? erm \ varies at different loca- 
tions of the compressed shells, e.g., more on the northeast 
(upper left) side and less in the southeast (lower left) di- 
rection for Run B. Within the shells surrounding the 
CR bubbles, the orientation of magnetic fields generally 
lies preferentially along the bubble surfaces. The align- 
ment and amplification is more pronounced for Run B 
(larger initial coherence length) than Run G, which ini- 
tially has a more randomly oriented field on small scales 
and the components perpendicular to the shock fronts 
are not draped effectively (e.g., the region at x ~ 4, 
—4 < z < —2 in the zoom-in image in Figure [3]). As we 
will show later, the magnetic field topology at the bub- 
ble edges has a critical impact on the appearance of the 
simulated bubbles. 

Inside the bubbles, the magnetic field is more uniform 
for Run B (Ib = 9 kpc), with an average strength of 
~ 0.1 //G, whereas for Run G (Ib = 1 kpc), the mag- 
netic field has a magnitude of ~ 0.1 — 1 /iG and is more 
filamentary. The enhanced field strength on the jet axis 
close to the GC comes from magnetic field injected by 
the AGN jets. The linear structure within the bubbles 
extending radially from the GC arises because of fields 
that are dragged along with the outwardly moving gas. 
We note that despite these enhancement, the magnetic 
field is dynamically unimportant for the expansion of the 
shocks and the CR bubbles, which is dominated by ther- 
mal pressure of the shocked gas and the CR pressure. In 
other words, the magnetic pressure is small compared to 
the total (gas plus CR) pressure, i.e., j3 > 1. We find 
that except in the magnetized jets and some filaments 
where plasma j3 ~ 1 — 5, the magnetic field is largely dy- 
namically negligible either in the CR bubbles (f3 > 50) 
or within the compressed shells (/3 ~ 10 — 50). Conse- 
quently, the CR distribution at t = ^Fermi for these two 
runs looks almost identical to the purely hydrodynamic 
case Run A (i.e., top panel of Figure [T]). 

The projected CR maps for simulations including CR 
diffusion are presented in Figure |4l Top and bottom pan- 
els are runs with diffusion coefficient 4 x 10 28 cm 2 s _1 
and 1.2 x 10 29 cm 2 s _1 , respectively. The panels from left 
to right show the cases of isotropic diffusion, anisotropic 
diffusion with Ib — 9 kpc, and anisotropic diffusion with 
Ib = 1 kpc, respectively. The isotropic cases with Ib = 1 
kpc (Run H and Run J) are omitted here because the re- 
sults are essentially identical to the corresponding Ib = 9 
kpc runs, because isotropic diffusion does not depend on 
the magnetic field geometry. Compared to the fiducial 
run without CR diffusion (Run A; Figure [2]), the cos- 
mic rays in the isotropic diffusion cases (Run C and Run 
E) have a smoother distribution because they are trans- 
ported by CR advection and additional diffusion. Be- 
cause of the dilution of CR energy density of the bubbles 
by diffusion, it takes a slightly longer time for the CR 
bubbles of e cr > 10~ n erg cm~ 3 to reach the observed 
latitude, i.e., ^Fermi = 1-22 Myr. Along the jet axis, the 
expansion of the CR bubbles is dominated by CR advec- 
tion (recall that the gas velocities is ~ 6000—7000 km s" 1 
in the vertical direction), while in the lateral direction, 
the relative contribution of diffusion is greater with re- 
spect to advection (the gas velocities in the lateral direc- 



tion is ~ 2500 km s _1 ). Consequently, the addition of CR 
diffusion has a more prominent effect laterally, yielding 
much fatter CR bubbles than the fiducial case. We also 
note that the diffusion transports cosmic rays to a larger 
distance than one would have estimated based on a sim- 
ple dimensional relation, I ~ y/nt Fermi, because of the 
large gradient of the CR energy density on the interface 
between the bubbles and the ambient gas. Therefore, the 
effect of CR diffusion is in fact quite substantial, and is 
easily visible in the figures. 

The most noticeable feature in the two cases with 
isotropic CR diffusion (Run C and Run E; left panels 
in Figure 2]) is the smooth transition of CR distribution 
at the edges, with CR energy density gradually declining 
outward. For k; so = 1.2 x 10 29 cm 2 s _1 , the smooth- 
ness due to diffusion can be easily seen, though may be 
exaggerated because such values of isotropic diffusion co- 
efficients may be hi gher than observationally constrained 
(|Strong et al.ll2007l) . However, even for an isotropic diffu- 
sion coefficient of 4x 10 28 cm 2 s _1 , as typically adopted in 
CR propagation m odels in the Galaxy (e.g., GALPROP; 
IStrong eta l. 2009), the resulting CR distribution still has 
a smooth edge in contradiction to the sharp edges seen in 
the observed gamma-ray map. The observed sharpness 
of bubble edges thus implies suppression of CR diffusion 
across the bubble surface. 

The maps of projected CR energy density for runs with 
anisotropic diffusion are shown in the middle and right 
columns in Figure [4] In contrast to the cases of isotropic 
diffusion, the edges of the CR bubbles are much sharper, 
similar to the fiducial run where CR diffusion is not in- 
cluded (except perhaps Run K). The sharp edges are 
attributed to the combined effects of magnetic draping 
and anisotropic CR diffusion. Because of the small gyro- 
radius of relativistic particles, CR diffusion is anisotropic, 
i.e., cosmic rays may diffuse primarily along magnetic 
field lines with strongly suppressed cross-field diffusion. 
During the formation of the CR bubbles, since the mag- 
netic field lines tend to align with the bubble surface due 
to magnetic draping (Figure[3]), CR diffusion occurs pref- 
erentially in the direction tangential to the bubble surface 
but is inhibited across the boundary. This suppression 
of CR diffusion across the bubble surface, which is a con- 
sequence of our MHD simulations that self-consistently 
include magnetic field and anisotropic CR diffusion, pro- 
vides a natural explanation to the sharp edges of the 
observed Fermi bubbles. 

We note that the presence of the magnetic draping 
effect is manifested by comparing the CR distribution of 
Run C (n iso = 4 x 10 28 cm 2 s" 1 ) and Run K (k\\ = 1.2 x 
10 29 cm 2 s^ 1 ). Since Run K has a tangled magnetic field 
on small scales, if the field is not aligned with the bubble 
surface but is completely random as it is initially, the 
result should correspond to an effective isotropic diffusion 
coefficient of ~ 1/3 of the parallel diffusion coefficient, 
i.e., Ki SO = 4 x 10 28 cm 2 s _1 as in Run C. The fact that 
the bubble edges are relatively sharper in Run K than 
in Run C indicates that, even though the magnetic field 
with a small coherence length may not be draped very 
effectively at some locations of the bubble surface, the 
magnetic draping effect does occur in a statistical sense 
(see the bottom panel of Figure [3]). 

Among the simulations with anisotropic diffusion (mid- 
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Fig. 4. — Maps of projected CR energy density in logarithmic scale in cgs units at t = tFermi for simulations including CR diffusion. Top 
and bottom rows show the runs with diffusion coefficient 4 X 10 cm 2 s" 1 (Run C, D, and I) and 1.2 X 10 29 cm 2 s" 1 (Run E, F, K), 
respectively. Columns from left to right are the cases of isotropic diffusion, anisotropic diffusion with l& = 9 kpc, and anisotropic diffusion 
with Zjri — 1 kpc, respectively. The solid and dashed lines are the observed northern and southern Fermi bubbles, respectively. Including 
isotropic diffusion results in a smoother CR distribution at the edges, whereas runs with anisotropic diffusion produce sharper edges like 
the observed Fermi bubbles. 



die and right columns in Figure U) , there are a few dif- 
ferences depending on the diffusion coefficients and the 
structure of the underlying GMF. If the diffusion coeffi- 
cient along the mean field is km = 4 X 10 28 cm 2 s _1 , the 
coherence length of the GMF has a minor effect on the 
appearance of the CR bubbles. On the other hand, if 
kii = 1.2 x 10 29 cm 2 s~ x , the (possibly) largest diffusion 
coefficient in the anisotropic diffusion case permitted by 
current observational limits (see § I2.1[) . the structure of 
the GMF may show observable features at the bubble 
edges. For example, in Run K (Ib = 1 kpc), where the 
magnetic field at the bubble surface is more randomly 
oriented (bottom right panel in Figure [3]), small-scale 
'wiggles' show up at the bubble edges due to fast diffu- 
sion along the field lines. The observed Fermi bubbles 
do not appear to have uneven surface at small scales. 
Therefore, by comparing our simulated CR image with 
the observation, we find that it is unlikely for our Galaxy 
to have fast diffusion (k|| ~ 1.2 x 10 29 cm 2 s _1 ) in a GMF 
that has a small coherence length (2b ~ 1 kpc). As for 
Run F (Ib = 9 kpc), the simulated northern CR bub- 
ble has a 'tail' near the Galactic plane toward the east 
(I ~ 10° - 20°, b ~ 0° - 5°), which reflects the orientation 
of the underlying magnetic field at this location (Figure 
[21 top right panel) . Although our initial magnetic field is 
likely not the same as the actual GMF, our simulation re- 
sults imply that the lack of absolute symmetry and even 
surface of the observed Fermi bubbles could be related to 
the structure of the underlying GMF. For instance, the 
slight extension or bend of the northern Fermi bubble 

at I 20°, b ~ 30° - 40° could possibly be the result 

of magnetic field that is oriented in the east-west (hori- 
zontal) direction and that was not effectively draped (see 
§ 14.21 for alternative explanations). Future observations 
of the GMF at the rims of the Fermi bubbles will help 



to verify whether this could be the case. 

For the Ib = 9 kpc cases shown in the middle col- 
umn of Figure |4l because the cosmic rays diffuse mainly 
along the field lines that are draped around the bub- 
ble surface, even a diffusion coefficient as large as km = 
1.2 x 10 29 cm 2 s" 1 is able to produce as sharp bubble 
edges as observed, just like the run with a smaller diffu- 
sion coefficient km = 4 x 10 28 cm 2 s _1 . This means that, 
if the GMF has a coherence length larger than the size 
of the Fermi bubbles, the CR diffusion coefficient par- 
allel to the mean field cannot be well constrained solely 
by the sharpness of the bubble edges. This is similar to 
fossil radio bubbles observed in galaxy clusters, where 
the cosmic rays may be confined within the bubbles due 
to magnetic drapi ng even when the para llel diffusivity 
coefficient is large (Ruszkowski ct al. 2008|). 

4. DISCUSSION 

4.1. Parameter Dependencies and Constraints 

In this section we comment on the parameter choices 
and dependencies. In the scenario where the Fermi 
bubbles are inflated b y jets from the central AGN, 
IGuo fc Mathewsl poll has performed a detailed survey 
of jet parameters and derived plausible parameters by 
matching the shape of the observed bubbles. We found 
the same parameter dependencies of bubble morphol- 
ogy as theirs - 'thinner' (vertically-elongated) bubbles 
are produced by jets with a larger density contrast n, a 
smaller total energy density in the jets ej + ej cr (which 
dominates the lateral expansion of the bubbles; also note 
that B 2 /8tt is negligible), a larger jet velocity Wj ct , a 
smaller jet width -R, e t, or a longer je t duration tj et . As 
discussed in IGuo &: Mathewsl (|201lD , these parameters 
are essentially degenerate in reproducing the observed 
bubble shape. Moreover, a number of jet parameters 
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derived from their simulations scale with the assumed 
normalization of the initial gas density profile (n e o), in- 
cluding |0j, ej, ej cr , and thus the total jet power Pj et and 
injected energy E- ]et . 

In order to reduce some of the parameter degenera- 
cies and improve the estimations of jet parameters, in 
this study we introduce several additional observational 
constraints. Firstly, for the hot gaseous halo we choose 
the initial temperature To and the normalization con- 
stant n fi o to matc h the observed gas density profile of 
iMiller fc Bregmanl ([20121 ) (§ This is particularly 

important because the supersonically expanding bubbles 
are partially confined by the pressure of the ambient gas, 
and thus the final shape of the bubbles depends on the 
chosen initial gas profiles. However, we note that there 
are still large uncertainties in the current co nstraints on 
the observed profiles of the hot gaseous halo (jFang et al.1 
120061 IMiller fc Bregmanl [2012] ) . Moreover, these mea- 
surements only provide information about the gas profile 
as observed today (especially after it is already modified 
by the AGN jet event in the inner few kpc), rather than 
the initial condition about ~ 1 Myr ago when the AGN 
jets were first launched. Therefore, though our simulated 
gas density profile at iFcrmi is within the observed limits, 
there are still uncertainties in the assumed initial hot gas 
distribution in the vicinity of the SMBH, and hence also 
in the estimated jet power and total energy. 

Secondly, observations of X-ray absorption line ratios 
of OVIII to OVII put a constraint on the temperature of 
the shocked gas inside the F ermi bubbles to be T > 10 7 
K (jMiller fc Bregmanl [20121) . We find that slower and 
wider jets are required to fulfill the constraints from the 
observed line ratios (§ I2.3p . implying that the relativis- 
tic jets from the AGN may have slowed down a little 
by interacting with the interstellar me dium during prop- 
agati on from pc to kpc scales (e.g., [Middclbcr g et al.l 
2004). Since these modifications would make the bub- 
bles much 'fatter', we had to either increase pj, ij ct , or 
decrease ej +ej cr in order to reproduce the observed bub- 
ble shape. However, p } and ij C t cannot be arbitrarily 
large, otherwise the electron number density at ^Fermi 
inside the bubbles would be too high and would vio- 
late the limb-brightened X-ray obse rvations by ROSAT 
([Snowden et al.lll997t Isii et al.ll2010D . Therefore, the to- 
tal energy density in the jets, ej + ej cr , becomes the pri- 
mary variable used to fit the observed shape of the bub- 
bles. 

We find that these observations set very stringent con- 
straints on the allowed parameter space for the jets; for 
a given initial gas profile, varying any of the jet parame- 
ters could easily violate one of the observational require- 
ments. It is remarkable that after applying these con- 
straints, the permitted parameters are able to simultane- 
ously reproduce the many characteristics of the observed 
Fermi bubbles, including the short age (thus consis- 
tent with the constraint from the IC cooling timescale), 
the absence of large-scale instabilities, and the limb- 
brightened CR distribution needed for uniform gamma- 
ray surface brightness. The only remaining degeneracy 
is between the injected thermal energy density ej and 
CR energy density ej cr , as their sum affects the pressure 
contrast with respect to the ambient medium and helps 
drive the expansion of the bubbles (though dominated 
by the ram pressure of the jets). This degeneracy can 



in principle be broken by detailed comparisons with ob- 
servations of radio synchrotron or IC emissions from CR 
electrons. 

4.2. Bends of the Bubbles 

The properties of our simulated CR bubbles are in 
broad agreement with the observed ones. However, the 
observed Fermi bubbles are not completely symmetric 
about the GC. Both the northern and southern bubbles 
are slightly bent to the west (negative longitude) . In this 
section we discuss possible causes of this asymmetry. 

Up to now we have only considered the simplest case 
of bipolar jets that are normal to the Galactic plane. 
However, AGN jets, particularly on kpc scales, are 
not necessarily aligned with the rotational axis of the 
galactic disks, and such misalig nments are observed fre- 
quently in other ga l axies (e.g., Middclb erg et all 120041: 
Giroletti et~aTI 120051: [Gallimore et al.ll2006t Kharb et al l 
20101) . and also recen tly found in the Milky Way 
(ISu fc Finkbemenl20Tl . " This misalignment could ei- 
ther arise from an intrinsic effect, su ch as jet precessions 
(e.g.. iFalceta-Goncalves et"al~ll2010h . or could be due to 
jet interactions with the surrounding medium as the jet 
travels from pc to kpc scales (e.g., Fanaroff- Riley type 
I radio sources). Forming the Fermi bubbles by intrin- 
sically tilted jets is appealing and is well-motivated by 
the s lanted gamma-ray jets recently discovered near the 
GC (|Su fc Finkbeinerl l2012'). However, in this case, the 
two jets would point at opposite directions at a 180° an- 
gle. They would potentially form two bubbles, with the 
northern one bending to the west, and the southern one 
to the east, inconsistent with the Fermi bubbles that 
both bend to the west. This argument, together with 
the fact that the gamma-ray jets have a spectrum at < 1 
GeV different from the Fermi bubbles (|Su fc Finkbeinerl 
12012ft , suggest that the bubbles probably do not originate 
from the newly-discovered gamma-ray jets. 

Bending of radio jets in some cases can be ex- 
plained by external ram press ure as galaxies move 
through an intracluster g as (e.g.. iBegelman et al.l 119791 : 
iBalsara k, Normanl IT992T ). It is known that the Milky 
Way moves in the intragroup medium (IGM). If the ram 
pressure of the wind generated by such motion is suffi- 
cient, it can potentially 'blow' the bubbles toward the 
same direction. Assuming the Fermi bubbles are tilted 
by a degrees from the rotational axis of the Galaxy, the 
required ram pressure from the wind ([Burns "et^fl979l) 
is Pram,wind = -Pram jet tana, where Pram jot = Pjvf ct is the 
ram pressure of the jets. Taking a = 10° and our fidu- 
cial jet parameters, P ra m,jct = 1-1 X 10~ 7 dyne cm~ 2 , we 
obtain P ram ,wind ~ 2x 10~ 8 dyne cm~ 2 , corresponding 
to a wind of density p w i n d ~ 4 x 10~ 24 g cm -3 and veloc- 
ity Wwind ~ 750 km s" 1 . These values are unreasonably 
large for estimated prope r ties of the Local Gr oup (e.g., 
IMcConnachie et al.l [20071 ICox fc Loebl l2008h . More- 
over, such strong wind woul d have dominated ove r the 
gravitational restoring force ([McCarthy et al.ll2008| ) and 
stripped a significant amount of Milky Way gas away. 
Therefore, given reasonable values of the IGM, the re- 
sulted wind is unlikely to have sufficient impacts on the 
bubbles to the desired degree. 

The required ram pressure, instead of being external, 
could possibly come from the region in the SMBH vicin- 
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ity. If the SMBH moves relative to the hot gaseous halo, 
it can effectively experience a wind and the ram pressure 
over time might affect the direction of the jets. How- 
ever, little motion of the SMBH at the GC is detected 
based on proper motion measurements of the Sgr A* 
(jReid k Brunthalerll2005l ). 

Another possibility is the ram pressure from supernova 
(SN) explosions. It is well known that the Sgr A* is sur- 
rounded by nuclear star clusters, among which exists a 
population of massive young stars that f ormed ~ 6 Myr 
ago i n the inner 0.5 pc from the SMBH (Paumar d et al.l 
2006) . These massive stars typically have main-sequence 
lifetimes of a few Myr, and hence it is possible that 
one of them exploded as a SN during or sometime be- 
fore the active phase of the AGN jets. In such a case, 
assuming a compressed interstellar medium of density 
p ~ 10 -25 g cm -3 and a typical SN shock velocity of 
v ~ 10 4 km s , the ram pressure generated by the SN 
would be sufficient to bend the jets. One could more ac- 
curately calculate the ram pressure as a function of time, 
the ambient gas density and the initial en ergy of explo - 
sion from, for example, the Sedov solution (|Sedovl fl946). 
We found that for a typical SN energy of 10 51 erg and 
ambient density of 10 -25 g cm~ 3 , the duration within 
which the SN can effectively bend the jets by more than 
10° is ~ 10 4 yr, which is about one tenth of the duration 
of our jets. 

In order to see whether such SN event could cause the 
observed bends of the Fermi bubbles of approximately 
the right morphology, we did a numerical experiment 
in which the northern and southern CR jets are both 
tilted by 10° to the negative x-axis (positive Galactic 
longitude) from the rotational axis of the Galaxy for 
< t < 0.1 tjet, and then returned to the normal axis for 
0.1 tjet < t < tjet- The projected CR energy density at 
^Fermi = 1.4 Myr is shown in Figure [5] As can be seen, 
the 'left '-bent jets result in a fatter CR distribution in 
the east, causing asymmetries about the rotational axis. 
The simulated bubbles have a remarkable resemblance to 
the observed bubble morphology, considering the simple 
assumptions and models employed. This suggests that 
bent AGN jets acting even for a short period of time, 
possibly due to ram pressure from a SN explosion that 
occurred near the SMBH, could plausibly cause the slight 
bends of the observed Fermi bubbles. 

As mentioned in § 13.21 the asymmetry of the observed 
bubbles could also be a result of large CR diffusivity 
combined with magnetic fields orienting perpendicular 
to the bubble surface. However, a special field configura- 
tion would be required for both bubbles to bend toward 
the same direction. Improved observations of the field 
geometry near the bubble edges will help to verify this 
possibility. 

4.3. Comparison with ROSAT X-ray Map 
The ROSAT X-ray 1.5 keV map (|Snowden et al.l 



Il997t ) has revealed enhanced emission surrounding the 
northern Fermi bubble, which is likely produced by 
bremsstra hlung emission f rom shocked gas during bubble 
formation (|Su et alJl20To() . The observed X -ray emission 
is limb-brightened, suggesting that the bubbles are hot 
and underdense. As discussed in § 14.11 this provides an 
important constraint on the thermal content of the AGN 
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Fig. 5. — Projected map of CR energy density in logarithmic 
scale in cgs units for the case where both CR jets are tilted to the 
east by 10° (possibly due to SN ram pressure; see texts) from the 
rotational axis of the Galaxy for < t < 0.1 tj ct , then returned 



to the normal direction for 0.1 t 



jet 



< t < tj 



jet- 



jet' 

Solid and dotted 

lines show the edges of the observed northern and southern Fermi 
bubbles, respectively. 



jets, i.e., the jets cannot be too 'heavy', otherwise the 
bubbles would be too bright on the X-ray map as they 
would be filled with large amounts of thermal gas. In 
§ 13.11 we have shown that our simulated bubbles are in- 
deed underdense and hot (see bottom panel of Figure 
[T]). However, the projections from the 3D distribution 
onto the 2D maps may be nontrivial, affecting the X-ray 
intensity distribution and location of the shocks, for in- 
stance. In order to show general consistency with the 
observed X-ray images, we make a simulated X-ray map 
by projecting the bremsstrahlung emissivity computed 
from the density and temperature of the simulated gas. 

For the simulated X-ray map, the X-ray emissivity 
in an energy range 1.4 — 1.6 keV is calculated using 
the MEKAL model felewe et al.l l l985tlKaastra k Mewej 



Il993t iLiedahl et al.l 19*95[ ) implemented in the utility 
XSPEC 4Arnau!l996) 

assuming solar metallicity. Note 
that the observed X-ray emission is contributed by all 
the gas in the Milky Way h alo, which likely extends 
to a radius of ~ 250 k pc ([Blitz k Robishawi 120001: 
iGrcevich k Putmanl l2009f ). much bigger than our sim- 
ulation box. Therefore, we first compute the X-ray emis- 
sivity from the simulated gas within a radius of 25 kpc 
away from the GC. Then, beyond 25 kpc the gas is as- 
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sumed to be isothermal with T — 2 x 10 6 K and fol- 
lows o ut to a radius of 25 kpc t he observed density pro- 
file of I Miller fc Bregmanl (|2012[ ) with an electron num- 
ber density floor, n e , floor- Studies of ram pressure strip- 
ping of dwarf spheroidal galaxies orbiting the Milky Way 
give constraints to the value of ro e , floor to be around 
2.4 x 1CT 5 - 2.5 x 10~ 4 c m" 3 ([Blitz fc Robishawl l200"ot 
iGrcevich fe Putnianl |2009() . Here we choose n e , floor = 
8.0 x 10~ 5 cm~ 3 so that the OVII column density com- 
puted from this exte nded gaseous halo is cons istent with 
the observed values ([Miller fc Brc gman 2012). We show 
the results from the same run as in the previous section 
(with a short period of jet bending) as it produces bub- 
bles with morphology closest to that observed. 

The projected X-ray emissivity at 1.5 keV is shown in 
Figure H2 overplotted with the observed contours of the 
Fermi bubbles and arc features identified by ROSAT. A 
region more extended than the bubbles is X-ray bright 
compared to the ambient medium because the gas is com- 
pressed and heated by the shocks during the bubble ex- 
pansion. The projected location of the shock fronts is in 
excellent agreement with the observed outer northern arc 
feature embracing the bubbles. The projected map has a 
much smoother distribution than one would have naively 
expected from a 2D slice that shows a clear cavity within 
the shock compressed shells (as in the bottom panel of 
Figure [1]). This is because the lines of sight toward the 
bubble interior also pass through the closer and further 
ends of the shells, and thus the change in projected X- 
ray emissivity with sight lines moving toward the shell 
edges is rather continuous. The projected X-ray emission 
inside the bubbles is dimmer and gradually increases to- 
ward the edges of the shocked gas, consistent with the 
limb-brightened distribution on the observed X-ray map. 
We note that the approach we use to produce the X-ray 
map is rather simplistic; simulated X-ray observations 
including stochastic photon generation and instrumen- 
tal responses for different bands are required to make 
more detailed comparisons with the observational data. 
However, the general success in reproducing the observed 
X-ray features provides a strong support for the hypoth- 
esis that both the Fermi bubbles and the ROSAT X-ray 
features originate from the same episode of jet activity 
from the central SMBH. 

5. CONCLUSIONS 

The Fermi Gamma-ray Space Telescope has recently 
revealed two large bubbles extending ~ 50° above and 
below the GC, with a width of ~ 40° in longitude. 
The northern and southern bubbles are nearly symmet- 
ric about the Galactic plane, with only slight bends to 
the west. Their spectrum is hard and has no signifi- 
cant variations within the bubbles. The gamma-ray sur- 
face brightness is quite uniform with sharp edges at the 
boundaries. Besides, the gamma-ray bubbles are spa- 
tially correlated with features in the ROSAT X-ray map 
and the hard-spectrum WMAP haze. It is challenging 
to explain all these observed properties, for example, by 
SN shoc ks in the Galac tic disk, or by dark matter anni- 
hilation (iSu et al.M20loh . 

The symmetry about the GC of the two bubbles, the 
bilobular shape, and their similar hard spectrum strongly 
suggest the bubbles originated from an energetic event in 
the GC within the past few Myr. Forming the bubbles 
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Fig. 6. — Simulated X-ray map at 1.5 keV (see texts for defini- 
tion) for the same run as Figure[B] The solid and dotted lines show 
the surfaces of the observed northern and southern Fermi bubbles, 
respectively. The dash lines are the inner and outer northern arcs 
observed by the ROSAT X-ray satellite. 



and their observed shape by a past jet a ctivity of the cen- 
tral SMBH is shown to be plausible bv lGuo fc Mathews! 
(2011) using 2D hydrodynamic simulations. They also 
found that the sharpness of bubble edges requires sup- 
pression of CR diffusion across the bubble surface, which 
may be related to the interplay between CR diffusion and 
magnetic fields. However, their simulated bubbles reveal 
several discrepancies with the observations, including the 
rippled surface due to hydrodynamic instabilities, and 
the inferred centrally-brightened surface brightness. To 
this end, shear v iscosity is invoked in their second paper 
(jGuo et al.ll2011| ) to alleviate these problems. 

In this study, we investigate the jet scenario using a set 
of 3D MHD simulations that self-consistently include the 
effects of magnetic fields and anisotropic (field-aligned) 
CR diffusion, as well as dynamical interaction between 
the thermal gas and cosmic rays. The simulations are 
performed using a new CR module in the FLASH code, 
which we have implemented and tested (see Appendix 
[A|) . We summarize our findings as follows. 

1. The effect of projection of the 3D bubbles has a 
significant impact on the estimation of the bubble for- 
mation time. Because the widths of the bubbles occupy 
a non-negligible fraction of the distance from the Sun to 
the GC, for the 3D CR bubbles to project onto a Galac- 
tic latitude of 50°, their vertical dimension only needs to 
be ~ 6 kpc, instead of ~ 10 kpc as previously thought. 
The projection effect results in a much shorter bubble 
formation time (~ 1.2 Myr for our fiducial model) than 
in previous estimations. 

2. If the observed gamma-ray emission is produced by 
IC scattering of CR electrons by the ISRF, the IC cool- 
ing time of high-energy (~ 100 GeV) electrons gives a 
stringent upper limit to the bubble ages to be within a 
few Myr. This constraint is naturally satisfied by the 
'young' bubbles revealed by our 3D simulations. 

3. Because of the short ages of our simulated bubbles, 
there is no sufficient time for large-scale hydrodynamic 
instabilities to grow. This alone explains why the ob- 
served Fermi bubbles have a rather smooth surface (as 
opposed to rippled) , and there is no need to invoke other 
mechanisms (e.g., viscosity or magnetic fields) to explain 
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the suppression of the instabilities. 

4. Our jet parameters, which ar e de termined by various 
observational constraints (see § 14. 1[) . produce an edge- 
enhanced CR distribution that, when projected onto the 
plane of the sky, result in a roughly flat surface bright- 
ness distribution as a function of the Galactic longitude. 
The projected CR intensity increases with the Galac- 
tic latitude, consistent with the flat gamma-ray surface 
brightness as the ISRF decreases away from the Galactic 
plane. 

5. The sharp edges of the Fermi bubbles are repro- 
duced by self-consistently including anisotropic CR dif- 
fusion along magnetic field lines that drape around the 
bubble surface during the bubble expansion. 

6. The causes of the slight bends of the Fermi bubbles 
are discussed. Possible explanations include jet bending 
due to ram pressure from SN near the SMBH, and fast 
CR diffusion along magnetic field lines that initially lie 
perpendicular to the bubble surface (see § 14. 2[) . 

7. The projected X-ray bremsstrahlung emissivity of 
the shocked gas due to bubble expansion successfully re- 
produces the location and limb-brightened property of 
the X-ray features surrounding the Fermi bubbles ob- 
served with ROSAT. This provides evidence that the 
ROSAT X-ray features originated from the same AGN 
jet activity episode as the Fermi bubbles. 

The properties of our simulated bubbles are in broad 
agreement with the key features of the Fermi bub- 
bles, strengthening the case that the bubbles are cre- 
ated by a recent AGN jet activity at the GC. Such 
event is analogous to bipolar radio bubbles and X-ray 
cavities associated with the nuclei of massive galaxies 
(|McNamara fc Nulsenl [2007). Fermi cavities are most 
likely created by jets from the central SMBH, e.g., as 
those observed in distant galaxies (e.g., M87). Though 
the SMBH at the center of the Milky Way is currently in 
a quiescent state (see, however, a pair of gamm a-ray jets 
recently discovered by ISu fc Finkbeinerl I2012D . the de- 
tection of the Fe rmi bubbles, together with indications 
of pa st activity (jKovama et al.l 119961 : iRevnivtsev et al.1 
2004), suggest that the central SMBH may have regularly 
undergone cycles of jet activity. The energy input from 
these jets may have altered the properties of the mul- 
tiphase gas in the Galactic bulge and halo, suppressed 
the accretion process, and regulated the co-evolution of 
the SMBH and the stellar bulge in the past. Moreover, 
these AGN jets may be a s ignificant source of the cosmic 
rays in the Galactic halo (|Guo fc Mathews! 120111 ). The 
proximity to the SMBH in the GC offers a unique op- 
portunity to study these processes in detail, especially in 
gamma rays due to the limited sensitivity and resolution 
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of gamma-ray observations. 

Finally, we note that in our simulations we have only 
considered the ensemble of cosmic rays and neglected 
distinctions between different CR species and energies. 
However, many physical processes depend on the types of 
CR particles under consideration and their energy spec- 
tra, such as CR production and reacceleration in shocks, 
and energy losses due to synchrotron and IC emission. 
Although our current results suggest that the observed 
features of the Fermi bubbles are reproduced reasonably 
well without these processes, it is possible that they also 
contribute to modify the properties of the bubbles. In 
particular, the observed Fermi bubbles have a hard spec- 
trum with little spatial variation within and between the 
two bubbles. Though the bubble formation time is quite 
short, as we have shown, the cooling of CR spectra due 
to synchrotron and IC radiation may still be relevant, 
depending on the CR energies and Galactic latitudes. 
Therefore, to produce the flat spectrum as observed is 
nontrivial and requires detailed modeling of the CR spec- 
tra. Furthermore, as discussed in § I4.H our current model 
is degenerate with respect to the internal energy density 
and CR energy density in the jets. A comparison of the 
simulated spectra in the gamma-ray and/or radio bands 
with the observational data should allow to break this 
degeneracy. This technique may also potentially offer a 
unique way to constrain the contents of the AGN bub- 
bles and is likely to be superior to the constraints from 
AGN bubbles in galaxy clusters because of the proxim- 
ity of the GC and the availability of additional spatially 
resolved data from Fermi. 
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APPENDIX 

A. NUMERICAL IMPLEMENTATIONS AND TESTS FOR COSMIC RAYS 

In this study, we introduce a new module in the FLASH code for simulating cosmic rays as a second fluid. A brief 
summary of the code and the directionally unsplit staggered mesh (USM) solver, as well as the general assumptions 
and equations for simulating CR advection and diffusion, are already discussed in § 12. H In this appendix, we first 
describe the implementation details of how we incorporate cosmic rays into the USM solver, and then provide results 
of several numerical tests of the new CR module. 

The MHD equations including CR advection, dynamical coupling of thermal gas and cosmic rays, and anisotropic 
CR diffusion, are listed in Eq. 1-6 in § 12.11 The mass, momentum, tot al energy, and inducti on equations are solved 
by th e Piecewise-Parabolic Method fPPM; ICoIeila fc Woodwardl (fl98l ) in the USM solver (jLee fc Deandl2009t iLei 
2012). In the absence of CR diffusion, the only modification is that the pressure contains an additional contribution 
from the cosmic rays, i.e., ptot = Pth + Pct +Pb = {l — l)eth + (7cr — l)e cr + B' 2 /8ir, where eth and e cr and the internal 
energy densities, and 7 and 7 cr are the adiabatic indices of the thermal gas and cosmic rays, respectively. Accordingly, 
the sound speed of the combined fluid of gas and cosmic rays, which is used when solving the Riemann problem in the 
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PPM method as well as when computing the hydrodynamic timestep, is replaced by the effective sound speed (|Miniatil 
[2003), 



7Pth + 7crPcr 
P 



(Al) 



We note that when the left and right states are reconstructed at cell boundaries in the PPM method, it is crucial that 
the CR energy is also reconstructed in the same way as other gas variables (i.e., to the same order and through the 
same slope limiters) in order to avoid spurious oscillations around discontinuities (e.g., Figure ITT]) . 

The energy of cosmic rays is evolved from Eq. [5J in which the advection term is advanced using the mass scalars in 
FLASH that passively evolve with the thermal gas, and then the source term — p cr V • v is updated in the discretized 
form (in 2D) as 



n+l/2 



n+1/2 



J i+l/2.j u i-l/2,j 

Ax 



V 



n+1/2 n+l/2 
i.J + 1/2 ~ V i,]-l/2 



Ay 



(A2) 



where Aa; and Ay are the sizes of grid cells, v n t}!n ■ and v n+ }( 2 /r , are velocities at cell boundaries reconstructed by the 

y & 5 i±l/2,j i,j±l/2 J 

PPM method, and p% i j is the averaged CR pressure before and after advection. Extension to 3D is straightforward, 
and we skip it here for the sake of brevity. 

The anisotropic CR diffusion terms in the total energy and CR energy e quations are updated explicitly using fluxes 
comp uted from the conservative 'centered asymmetric differencing' scheme (Parrish & Stone 2005; Sharma fc Hammettj 
2007). The updated energy (for either the total or CR energy) in 2D is 



„n+l 



At 



r i+l/2,J 



zpn 



r i,j+l/2 



F, 



(A3) 



Aa; Ay 

The simulation timestep At is required to satisfy both the CFL criterion and the stability condition for CR diffusion, 

min[Aa; 2 , Ay 2 } 



At < 



2(k 



(A4) 



Using the asymmetric method, the x-flux (similar for the y- fluxes) at the cell face i + 1/2 at time n due to anisotropic 
CR diffusion (Eq. [5J) is given by 



bm = 



b x 



de cr 
dx 



- de cr 
'~dy~ 



, de CI 
dx ' 



u x,i+l/2,ji 
de CY ^cr,«+l,j &cx,i,j 

dx Ax 

by = (by,i,j-X/2 + byS+1 j-l/2 +^i/,ij + l/2 + ^jm+I J+l/2)/4, 
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(A6) 
(A7) 
(A8) 
(A9) 



and the transverse flux is s lopcd-limited in order to avoid negative CR energy densities when there is a large gradient 
(Sh arma fc Hammettll2007ft . 
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= L{L 
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de r . 
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where L is a slope limiter such as minmod, van Leer, or Monotonized Central (MC) limiter. We adopted the MC 
limitcr, defined as 



MC(a, b) = minmod 
minmod(a, b) 



a + b 



2minmod(a, b), 



min(a, b) if a, b > 0, 
max(a, b) if a, b < 0, 
if ab < 0. 



(All) 
(A12) 



A.l Test for CR Advection 

First we perform a simple 2D test where cosmic rays are advected along the diagonal direction. The simulation 
domain of size 1 x 1 is initially filled with a uniform thermal gas of density po and pressure po, moving with the 
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Fig. 7. — Test for CR advection. Left and right panels are the CR energy density at t = and t = 0.25 s. Adaptive mesh tracks the 
overdensity of cosmic rays without introducing spurious features. 



velocity (v x o, v y o). On top of the thermal gas we placed a Gaussian overdensity of cosmic rays, 

e C r = e cr , (l + cxp[-0.5(r/r ) 2 ]). (A13) 

The following parameters are adopted: po = 1, po = 1, e cr .o = 0.1, ro = 0.05, v x q = 1, and v y o = 1. Since the purpose 
of this test is only to verify the implementation of CR advection, the CR pressure is set to be much smaller than the 
thermal gas pressure in order to minimize the back reaction of cosmic rays on the gas. The simulation is performed 
with adaptive mesh refinement (AMR), with a minimum and maximum resolution of 1/256 and 1/2048, respectively, 
depending on the local second derivative of the CR energy density. Magnetic fields and CR diffusion are not included 
in this run. 

Figure [7] shows the CR energy density at t = (left) and t — 0.25 s (right), overplotted with the locations of the 
AMR blocks (each block contains 8x8 grid cells). The cosmic rays passively evolve with the thermal gas along the 
diagonal direction without any distortion. Moreover, this test demonstrates that the AMR correctly traces the CR 
overdensities and does not introduce any artifacts. 

A. 2 Tests for Magnetic Field Aligned CR Diffusion 

In order to verify the component of anisotropic CR diffusion in the new CR module, we put an overdensity of cosmic 
rays (Eq. IA13I) in two different magnetic field configurations: (1) diagonal field, and (2) loop field. In order to focus 
on the effect of CR diffusion, in these tests we turn off the calculation of hydrodynamic fluxes and only solve for the 
CR evolution due to diffusion along the field lines according to Eq. [6] CR back reaction on the magnetic field is also 
neglected. 

For the first 'diagonal' test, the initial CR distribution is the same as in the advection test (see § IA.1[) . except that 
the fluids are at rest in the beginning of the simulation, and the cosmic rays are allowed to diffuse along the magnetic 
field with diffusion coefficients «|| = 0.05, and kj_ = 0. Here we present results for a simulation box of 1 x 1 with 
a uniform grid of 256 x 256 cells. Figure [8] displays the evolution of CR energy density in a diagonal magnetic field 
(arrows). The cosmic rays diffuse in the diagonal direction along the field lines as expected. 

Next, we show the test results of anisotropic CR diffusion along a 'looped' magnetic field, which is more stringent 
because the field lines are inclined at all possible angles with respect to the Cartesian grid. For this test, we initialize 
a magnetic field with uniform amplitude Bo — 10~ 4 /iG in concentric circles around the center of the simulation box. 
This test is performed in 2D, and a virtual vector potential (B = V x A) is used to ensure V • B — 0. To generate a 
uniform looped field, the vector potential is chosen to be A — (A x , A y , A z ) — (0,0, i?o(5 — r)), where r is the radius 
from the center of the domain. The initial CR distribution is the same as in the diagonal test, but displaced by 0.1 
from the center of the simulation box. The CR distribution at three different epochs is shown in Figure [9l Due to 
anisotropic diffusion around the looped field (arrows), the cosmic rays start from a initially localized distribution and 
gradually diffuse following the field lines, and eventually equilibrate when the CR energy gradient in the azimuthal 
direction vanishes. This stringent loop test demonstrates the robustness of our implementation of CR anisotropic 
diffusion. 



19 













time 


= 




o.ooooc 


ps 




















time 






1. 




s 










1 .0 


? 


/ 


/ 


1 

A 


/ 


/ 




./ 


/ 


/ 


A 


/ 


/ 


a 


/ 


/ ■ 




1 .0 




/ 


y 


y 


/ 


/ 


1 

A 


.A 


/ 


A 


y 


/ 


A 


A 


A 


A - 




A 


/■ 


A 


A 


A 


A 


A 


A 


A 


A 


A 


A 


a 


a 


y 


/• 






y 


/ 


/■ 


y 


y 


/ 


y 


A 


A 


A 


A 


A 


A 


A 


A 


A ' 


0.3 


A 


A 


A 


A 


./ 


A 


A 


A 


A 


/ 


A 


A 


a 


a 


A 


' 




0,0 


> 


/ 


y 


y 




A 


A 


A 


■A 


A 


y 


A 


A 


A 


A 


A ' 




r 


/ 


/ 


A 


/ 


./ 


/ 


/ 


/ 


/ 


/ 


A 


/ 


/ 


f 


/ . 




/ 


/ 


/ 


/ 


/ 


/ 


/ 


/ 


/ 


/ 


A 


A 


A 


A 


/ 


/ . 




V 


/ 


A 


f 


/ 


/ 


A 


A 


A 


/ 


A 


A 


/ 


a 


f 


/ ■ 






-/ 


/ 


y 


y 


/ 


y 


A 


A 


A 






A 


A 


A 


A 


A - 




a 


A 


A 


a 


/ 


A 


A 






A 


A- 


A 


a 


a 




/• " 






-/ 




y 


y 


A 


A 


A 


A 










A 


A 


A 


A ~ 


0.6 


> 


A 


A 


a 


A 


A 


A 










A 


a 


a 




y 




0.6 


> 




y 


y 


y 


A 


A 










f 


A 


A 


A 


A^ 












A 


A 




d 




\ 


\ 


/ 


/ 


















/ 


/ 
























<* 


A 


A 


A 


/ 


/ 


s 


i 


I 






A 


a 


/ 


/ 


/ ■ 






/ 


/ 


y 


./ 


/ 












A 


A 


A 


/ 


A 


A - 


0.4 




A 


A 


a 


A 


a 












A 


/ 


a 


J* 






0.4 




y 


y 


y 


y 










A 


A 


A 


A 


A 


A 


A~ 




> 


A 


A 


a 


/ 


A 


A 






A 


A 


A 


a 


A 




/• ' 






> 


/ 


y 


y 


A 










A 


A 


A 


A 


A 


A 


A " 




** 


A 


A 


A 


/ 


A 


■J 


A 


/ 


/ 


A 


A 


/ 


/ 


/ 








/ 


/ 




/ 


/■ 


I 


I 


w 


/ 


A 


J 


A 


A 


A 


A 


A . 


0,2 




/ 


A 


A 


/ 


f 


A 


J 


A 


/ 


J* 


/ 


/ 


/ 


f 


/- 




0.2 


? 


/ 


y 


y 


/ 


A 


A 


A 


A 


A 


A 


A 


A 


A 


A 


A~ 






A 


A 


A 


A 


A 


a 


a 


A 


? 


A 


A 


a 


A 




/• ■ 






y 


y 


y 


y 


A 


A 


A 


y 


A 


y 


y 


A 


A 


A 


A 


A ~ 






A 


A 


a 


A 


a 


a 


a 


A 


A 


A 


A 


a 


A 




/ 






> 


y 


y 


y 


y 


A 


y 


y 


A 


A 


A 


A 


A 


A 


A 


A ' 


0.0 




/ 


A 


If 


A 


/ 


A 




A 


/I 




r 


,/ 




A 






0.0 




? 


A 


j? 


,/ 




A 




A 


, F\ 


A 


/ 


A 




A 


S 



; (cm) 



; (cm) 



FlG. 8. — Evolution of the CR energy density due to diffusion along a diagonal magnetic field (arrows). 
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Fig. 9. — Time series of the distribution of the CR energy density. The cosmic rays gradually diffuse along the looped magnetic field 
(arrows) . 
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FlG. 10. — Linear sound wave test for a composite of thermal gas and cosmic rays. The figure shows the perturbation of the CR pressure 
along the x-direction at y = 0.5 after one wave crossing time. Simulated data is plotted with the plus symbols; solid line shows the 
analytical solution. 
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A. 3 Linear Sound Wave Test for a Composite of Thermal Gas and Cosmic Rays 

Next, we perform additional quantitative tests of the mod ule. Here we present result s of a 2D test of a sound wave 
propagating in a composite of thermal gas and cosmic rays (jRasera fc Chan dran 2008). This test problem takes into 
account the coupling between the two fluids. The medium is perturbed by small fluctuations that satisfy the following 
relations 

8 -lJ-l, (A14) 
Po c s 

Sp Sv ,. ir , 
— =7— , (Al5 

PO C s 

=7cr — , (Al6) 

Pcr,0 C s 

where Sp, Sv, Sp, and Sp cr and perturbations, p , p , p cr are the initial unperturbed quantities, and the adiabatic 
wave speed is 

IPO + 7crPci~0 , Al7 s 
PO 

For our test, we initialize a sine wave moving in the a;-direction with wavelength 0.5 and amplitude Sv = 10~ 3 in 
a medium with unperturbe d quantiti es po — l, vq — 0, po — l, and p C r.o — I- The other perturbed variables are 
calculated according to Eq. IAl5llAl6l All variables have values independent of their y coordinates. The wave travels 
periodically in a simulation box of sizes l x l with 128 grid cells on a side. Figure [TOl shows the results after propagating 
the sound wave for one period. Our results (plus signs) show excellent agreement with the analytical solution (solid 
line). 

A. 4 Nonlinear Shock Tube Test for a Composite of Thermal Gas and Cosmic Rays 

The Sod shock tube problem (Sod 1978) is a standard test for the accuracy of computational fluid codes, in particular 
the Riemann solver. The original problem contains a polytropic gas separated into a high pressure and high density 
left state and a low pressure and low density state on the right side, which then evolves into a system characterized 
by a rarefaction wave, a contact discontinuity, and a shock. The analytical solution for the propagation of these 
characteristics can be derived using the Rankine-Hugoniot jump conditions. 

Similarly, for the case of a hybrid fl uid consisting of therm al gas and cosmic rays, an analytical solution of the shock 
tube problem has been derived by iPfrommer et all (j2006) . We perform this test in 2D in a simulation box of size 
1 x 1 on a uniform grid with 1024 cells on a side. The variables of the left state (0 < x < 0.5) are initialized to be 
Pl = 1, vl = 0, pl = 6.7 x 10 4 , and p cr ,L = 1-3 x 10 5 ; the right state (0.5 < x < 1) has values of pr = 0.2, Vr = 0, 
Pr — 2.4 x 10 2 , and p CI: R — 2.4 x 10 2 . The values on the cells are independent of their y coordinates. 

The results after running the simulation for t = 4.4 x 10 -4 are plotted in Figure [TT] Again, the simulated results 
(plus signs) match very well with the analytical prediction (solid line) . The transitions between characteristics are well 
located. The profile of the rarefaction wave is reproduced, and the contact discontinuity and the shock are resolved 
within few cells without spurious oscillations. Based on the above test results, we conclude that our implementation 
of the CR module in FLASH is successful. 
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Fig. 11. — Shock tube test for a hybrid fluid of thermal gas and cosmic rays. The figure shows profiles of the density (top left), velocity 
(top right), gas internal energy density (bottom left), and CR energy density (bottom right) sliced through y = 0.5. Plus signs are the 
simulated data, solid lines are the analytical solution, and dotted lines show the initial condition. 



